using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// DQC25C - Integration rules for the computation of Cauchy principal value integrals.
    /// </summary>
    [Serializable]
    public static class Dqc25cClass
    {
        /// <summary>
        /// Integration rules for the computation of Cauchy principal value integrals.
        /// </summary>
        /// <param name="f">The double precision function defining the integrand.</param>
        /// <param name="a">The left end point of the integration interval.</param>
        /// <param name="b">The right end point of the integration interval.</param>
        /// <param name="c">The parameter in the weight function.</param>
        /// <param name="abserr">The estimate of the modulus of the absolute error.</param>
        /// <param name="krul">Key which is decreased by 1 if the 15-point Gauss-Kronrod scheme is used.</param>
        /// <param name="neval">The number of function evaluations.</param>
        /// <returns>The result of the Cauchy principal value integral.</returns>
        public static double Dqc25c(HardRealFunction f, double a, double b, double c, ref double abserr, ref int krul,
                                    ref int neval)
        {
            double[] x = {
                             0.99144486137381041114, 0.96592582628906828675, 0.92387953251128675613,
                             0.86602540378443864676
                             , 0.79335334029123516458, 0.70710678118654752440, 0.60876142900872063942,
                             0.50000000000000000000, 0.38268343236508977173, 0.25881904510252076235,
                             0.13052619222005159155
                         };

            double ak22;
            double amom0;
            double amom1;
            double amom2;
            double cc;
            double centr;
            double[] cheb12 = new double[13];
            double[] cheb24 = new double[25];
            double[] fval = new double[25];
            double hlgth;
            double p2 = 0;
            double p3 = 0;
            double p4 = 0;
            double resabs = 0;
            double resasc = 0;
            double res12;
            double res24;
            double u;
            double result;
            int i;
            int isym;
            int k;
            int kp = 0;

            cc = (2.0 * c - b - a) / (b - a);
            if (Math.Abs(cc) < 1.1)
            {
                goto _10;
            }

            //  Apply the 15-point Gauss-Kronrod scheme.    
            (krul)--;
            result = Dqk15wClass.G_K15W(f, DqwgtClass.Dqwgtc, c, p2, p3, p4, kp, a, b, ref abserr, ref resabs,
                                        ref resasc);
            neval = 15;
            if (resasc == abserr)
            {
                (krul)++;
            }
            goto _50;

            //  Use the generalized Clenshaw-Curtis method. 
            _10:
            hlgth = 0.5 * (b - a);
            centr = 0.5 * (b + a);
            neval = 25;
            fval[0] = 0.5 * f.SolveAt(hlgth + centr);
            fval[12] = f.SolveAt(centr);
            fval[24] = 0.5 * f.SolveAt(centr - hlgth);
            for (i = 1; i < 12; i++)
            {
                u = hlgth * x[i - 1];
                isym = 24 - i;
                fval[i] = f.SolveAt(u + centr);
                fval[isym] = f.SolveAt(centr - u);
            }

            //  Compute the Chebyshev series expansion. 
            DqchebClass.Dqcheb(x, fval, cheb12, cheb24);

            //  The modified Chebyshev moments are computed by forward
            // *  recursion, using amom0 and amom1 as starting values.
            // 
            amom0 = Math.Log(Math.Abs((1.0 - cc) / (1.0 + cc)));
            amom1 = 2.0 + cc * amom0;
            res12 = cheb12[0] * amom0 + cheb12[1] * amom1;
            res24 = cheb24[0] * amom0 + cheb24[1] * amom1;
            for (k = 2; k < 13; k++)
            {
                amom2 = 2.0 * cc * amom1 - amom0;
                ak22 = (k - 1) * (k - 1);
                if ((k / 2) * 2 != k)
                {
                    amom2 -= (4.0 / (ak22 - 1.0));
                }
                res12 += (cheb12[k] * amom2);
                res24 += (cheb24[k] * amom2);
                amom0 = amom1;
                amom1 = amom2;
            }
            for (k = 13; k < 25; k++)
            {
                amom2 = 2.0 * cc * amom1 - amom0;
                ak22 = (k - 1) * (k - 1);
                if ((k / 2) * 2 != k)
                {
                    amom2 -= (4.0 / (ak22 - 1.0));
                }
                res24 += (cheb24[k] * amom2);
                amom0 = amom1;
                amom1 = amom2;
            }
            result = res24;
            abserr = Math.Abs(res24 - res12);
            _50:
            return result;
        }
    }
}